#include "numericalintegrator.h"
#include <QString>

double NumericalIntegrator::calculateArea(const std::vector<double>& x,
                                          const std::vector<double>& y,
                                          IntegrationMethod method)
{
    QString errorMsg;
    if (!validateData(x, y, errorMsg)) {
        throw std::runtime_error(errorMsg.toStdString());
    }

    switch (method) {
    case TrapezoidalRule:
        return trapezoidalRule(x, y);
    case SimpsonsRule:
        return simpsonsRule(x, y);
    case BoolesRule:
        return boolesRule(x, y);
    default:
        throw std::runtime_error("Unknown integration method");
    }
}

bool NumericalIntegrator::validateData(const std::vector<double>& x,
                                       const std::vector<double>& y,
                                       QString& errorMessage)
{
    if (x.size() != y.size()) {
        errorMessage = "X and Y arrays must have the same size";
        return false;
    }

    if (x.size() < 2) {
        errorMessage = "At least 2 data points are required";
        return false;
    }

    // 检查x是否单调递增
    for (size_t i = 1; i < x.size(); ++i) {
        if (x[i] <= x[i-1]) {
            errorMessage = "X values must be strictly increasing";
            return false;
        }
    }

    return true;
}

bool NumericalIntegrator::isEquallySpaced(const std::vector<double>& x, double tolerance)
{
    if (x.size() < 2) return false;

    double dx = x[1] - x[0];
    for (size_t i = 2; i < x.size(); ++i) {
        if (std::abs((x[i] - x[i-1]) - dx) > tolerance) {
            return false;
        }
    }
    return true;
}

double NumericalIntegrator::trapezoidalRule(const std::vector<double>& x, const std::vector<double>& y)
{
    double area = 0.0;
    for (size_t i = 1; i < x.size(); ++i) {
        double dx = x[i] - x[i-1];
        double avgHeight = (y[i] + y[i-1]) / 2.0;
        area += avgHeight * dx;
    }
    return area;
}

double NumericalIntegrator::simpsonsRule(const std::vector<double>& x, const std::vector<double>& y)
{
    if (x.size() < 3 || x.size() % 2 == 0) {
        throw std::runtime_error("Simpson's rule requires odd number of points (>=3)");
    }

    if (!isEquallySpaced(x)) {
        throw std::runtime_error("Simpson's rule requires equally spaced x values");
    }

    double dx = x[1] - x[0];
    double sum = y[0] + y[y.size()-1];

    for (size_t i = 1; i < y.size()-1; ++i) {
        if (i % 2 == 1) {
            sum += 4.0 * y[i]; // 奇数项
        } else {
            sum += 2.0 * y[i]; // 偶数项
        }
    }

    return sum * dx / 3.0;
}

double NumericalIntegrator::boolesRule(const std::vector<double>& x, const std::vector<double>& y)
{
    if (x.size() < 5 || (x.size() - 1) % 4 != 0) {
        throw std::runtime_error("Boole's rule requires 4n+1 points (n>=1)");
    }

    if (!isEquallySpaced(x)) {
        throw std::runtime_error("Boole's rule requires equally spaced x values");
    }

    double dx = x[1] - x[0];
    double sum = 0.0;

    for (size_t i = 0; i < x.size() - 1; i += 4) {
        if (i + 4 < x.size()) {
            sum += 7 * y[i] + 32 * y[i+1] + 12 * y[i+2] + 32 * y[i+3] + 7 * y[i+4];
        }
    }

    return 2.0 * dx * sum / 45.0;
}

double NumericalIntegrator::linearInterpolation(const std::vector<double>& x,
                                                const std::vector<double>& y,
                                                double targetX)
{
    if (targetX <= x.front()) return y.front();
    if (targetX >= x.back()) return y.back();

    for (size_t i = 1; i < x.size(); ++i) {
        if (targetX <= x[i]) {
            double t = (targetX - x[i-1]) / (x[i] - x[i-1]);
            return y[i-1] + t * (y[i] - y[i-1]);
        }
    }

    return y.back();
}
